% Axisymmetric KZK simulator for nonlinear sound beams from piston
% transducer (according to Lee, Y.S., et al., 'Time-domain modeling of
% pulsed finite-amplitude sound beams', J Acoust Soc Am, vol. 97, no. 2,
% 1995).
%
% Temperature dependend values are set for water at 22.5 degrees Celcius.
%
% author:      Marco M. Voormolen
% version:     1.0
% first draft: 30 March 2007
% last update: 6 May 2010


clear all
close all
fclose('all');


% Input parameters

  % Transducer
Radius=0.75*25.4e-3/2; % transducer radius in m

  % Excitation
p_0=10500*3.2; % the source pressure amplitude in Pa
f_0=1e6; % center frequency of the excitation pulse in Hz
FracBW=30; % fractional bandwidth of the excitation pulse in %

  % Calculation (spatial and temporal discretization)
N_Sample=80; % relative sample frequency in relation to the excitation frequency (f_0)
N_Apert=50; % number of points across the radius of the transducer
dSigmaIBFD=0.2e-3; % step size in the propagation direction for the dimensionless axial coordinate used for the IBFD method
dSigmaCNFD=2e-3; % step size in the propagation direction for the dimensionless axial coordinate used for the CNFD method
z_End=90e-3; % in m
% Obtain lateral profiles by setting the 'z_End' at the desired depth.

  % Results
PlotHarBln=1; % plot the first three harmonic components or only the fundamental
FiltOrd=5; % filter order of the Butterworth filter to obtain the amplitude of the envelope for the various components
FiltBW=[-0.5 0.5]; % relative bandwith of the Butterworth filter for the various components ('([<lower bound> <upper bound>] + i)*f_0')


% Calculated parameters
Omega_0=2*pi*f_0;
c_0=1489.69;
z_0=(Omega_0*Radius^2)/(2*c_0);
Alpha_0=23.20e-3*(f_0/1e6)^2;
A_Absorp=Alpha_0*z_0;
Rho_0=997.67;
Beta=1 + 0.5*5.03;
z_Bar=(Rho_0*c_0^3)/(Beta*Omega_0*p_0);
N_Nonlin=z_0/z_Bar;
N_Cycle=100*(4*sqrt(log(2)))/(pi*FracBW);
N_tMin=2*N_Sample*N_Cycle;
N_tMax=4*N_Sample*N_Cycle;
f_Sample=N_Sample*f_0;
t=-N_tMin/f_Sample:1/f_Sample:N_tMax/f_Sample;
N_t=length(t);
dTau=(2*pi)/N_Sample;
Tau=-N_tMin*dTau:dTau:N_tMax*dTau;
dRho=1/(N_Apert + 0.5);
N_r=3*N_Apert;
i=1;
SigmaMat(i)=0;
while SigmaMat(i)<z_End/z_0
  if SigmaMat(i)<=0.1
    dSigmaMat(i)=(1 + SigmaMat(i))^2*dSigmaIBFD;
  else    
    dSigmaMat(i)=(1 + SigmaMat(i))^2*dSigmaCNFD;
  end    
  i=i + 1;  
  SigmaMat(i)=SigmaMat(i - 1) + dSigmaMat(i - 1);  
end
N_z=length(SigmaMat);
dr_0=dRho*Radius;
r_0=[0:dr_0:(N_r - 1)*dr_0];
P=zeros(N_t, N_r);
for j=1:N_Apert + 1
  Excitation=exp(-(((Tau + ((j - 1)*dRho)^2)/(pi*N_Cycle)).^2)).*-sin(Tau + ((j - 1)*dRho)^2);
  P(:, j)=Excitation';
end
P([1 end], :)=0;
P(:, end)=0;
P_Source=P;
P_Source(1, end)=0;


% Propagate the source field
p_Axial=zeros(N_t, N_z);
z_Axial=zeros(1, N_z);
h=waitbar(0, ['Progress (z=' num2str(z_Axial(1)*1e3, '%.2f') 'mm)']);
for i=1:N_z - 1

  Sigma=SigmaMat(i);

  p_Axial(:, i)=squeeze(p_0*P(:, 1)/(1 + Sigma));
  z_Axial(i)=Sigma*z_0;
    
  dSigma=dSigmaMat(i);  
  if Sigma<=0.1
    DiffracIBFD(int32(N_t), int32(N_r), Sigma, dSigma, dTau, dRho, P)
    AbsorpIBFD(int32(N_t), int32(N_r), dSigma, dTau, A_Absorp, P)
  else  
    DiffracCNFD(int32(N_t), int32(N_r), Sigma, dSigma, dTau, dRho, P)
    AbsorpCNFD(int32(N_t), int32(N_r), dSigma, dTau, A_Absorp, P)
  end
  tic
  P=Nonlin_solution(Sigma, dSigma, dTau, N_Nonlin, N_t, N_r, P);
  waitbar(i/(N_z - 1), h, ['Progress (z=' num2str(z_Axial(i)*1e3, '%.2f') 'mm)'])
  
end
p_Axial(:, N_z)=squeeze(p_0*P(:, 1)/(1 + SigmaMat(N_z)));
z_Axial(N_z)=SigmaMat(N_z)*z_0;
close(h)


% Plot results

  % Analytical
t_Ana=-2*N_Cycle/f_0:1/f_Sample:6*N_Cycle/f_0;
p_Ana=zeros(length(t_Ana), N_z);
p_Front=exp(-(((Omega_0*t_Ana)/(pi*N_Cycle)).^2)).*-sin(Omega_0*t_Ana);
for i=1:N_z
  t_Edge=t_Ana - (sqrt(z_Axial(i)^2 + Radius^2) - z_Axial(i))/c_0;
  p_Edge=exp(-(((Omega_0*t_Edge)/(pi*N_Cycle)).^2)).*-sin(Omega_0*t_Edge);
  p_Ana(:, i)=p_0*(p_Front - p_Edge)';
end
figure(1)
h_1(1)=plot(z_Axial*1e3, 20*log10(max(abs(hilbert(p_Ana)))), 'b-', 'LineWidth', 2);
figure(2)
N_FFT=max([2^ceil(log2(length(t_Ana))) 2^ceil(log2(N_t))]); 
h_2(1)=plot(z_Axial*1e3, 20*log10(max(abs(fft(cat(1, p_Ana, zeros(N_FFT - length(t_Ana), N_z)))/f_Sample))), 'b-', 'LineWidth', 2);

  % KZK  
figure(1) % Axial profile (maximum amplitude of the envelope)
hold on
if PlotHarBln
  P_Filt=zeros(3, N_z);
  for i=1:3
    [B, A]=butter(FiltOrd, ((FiltBW + i)*f_0)/(f_Sample/2));
    DummyMat=filter(B, A, p_Axial);
    P_Filt(i, :)=20*log10(max(abs(hilbert(DummyMat))));
    h_1(2)=plot(z_Axial*1e3, P_Filt(i, :), 'r-', 'LineWidth', 1);
  end
  set(gca, 'YLim', [-60 0] + ceil(max(P_Filt(1, :))/5)*5)
else  
  h_1(2)=plot(z_Axial*1e3, 20*log10(max(abs(hilbert(p_Axial)))), 'r-', 'LineWidth', 1);
end
box on
xlabel('\bf{z} \rm{[mm]}')
ylabel('\bf{p} \rm{[dBPa]}')
title('maximum amplitude of the envelope')
legend(h_1, {'analytical'; 'KZK'}, 'Location', 'SouthEast')

figure(2) % Axial profile (maximum amplitude of the spectrum)
hold on
if PlotHarBln
  df=f_Sample/N_FFT;
  f=0:df:floor(N_FFT/2)*df;
  p_Spec=abs(fft(cat(1, p_Axial, zeros(N_FFT - N_t, N_z)))/f_Sample);
  P_FFT=zeros(3, N_z);
  for i=1:3
    [A, n_Lower]=max(f>(i - 0.5)*f_0);
    [A, n_Upper]=max(f>=(i + 0.5)*f_0);
    P_FFT(i, :)=20*log10(max(p_Spec(n_Lower:n_Upper - 1, :)));
    h_2(2)=plot(z_Axial*1e3, P_FFT(i, :), 'r-', 'LineWidth', 1);
  end
  set(gca, 'YLim', [-60 0] + ceil(max(P_FFT(1, :))/5)*5)
else  
  h_2(2)=plot(z_Axial*1e3, 20*log10(max(abs(fft(cat(1, p_Axial, zeros(N_FFT - N_t, N_z))))/f_Sample)), 'r-', 'LineWidth', 1);
end
box on
xlabel('\bf{z} \rm{[mm]}')
ylabel('\bf{p} \rm{[dB]}')
title('maximum amplitude of the spectrum')
legend(h_2, {'analytical'; 'KZK'}, 'Location', 'SouthEast')

figure(3) % Envelope of the last 'P(t, r)' plane
N_tResamp=max([1 floor(N_t/200)]);
N_rResamp=max([1 floor(N_r/100)]);
P_Dummy=P(1:floor(N_t/200):end, 1:floor(N_r/100):end);
dr=dRho*Radius*(1 + SigmaMat(N_z));
N_rDummy=size(P_Dummy, 2);
x=[-(N_rDummy - 1):N_rDummy - 1]*N_rResamp*dr*1e3;
y=t(1:N_tResamp:end)*1e6;
imagesc(x, y, abs(hilbert(cat(2, flipdim(P_Dummy(:, 2:end), 2), P_Dummy))))
axis tight
box on
xlabel('\bf{r} \rm{[mm]}')
ylabel('\bf{t'' - (r/R)^2/(1 + \sigma)/\omega_0} \rm{[µs]}')
title('envelope')
x_Lim=get(gca, 'XLim');
figure(3)
set(gca, 'XLim', x_Lim)
  
figure(4) % Lateral profile (maximum amplitude of the envelope)
hold on
if PlotHarBln
  P_Filt=zeros(3, N_r);
  for i=1:3
    [B, A]=butter(FiltOrd, ((FiltBW + i)*f_0)/(f_Sample/2));
    DummyMat=filter(B, A, P*p_0/(1 + Sigma));
    P_Filt(i, :)=20*log10(max(abs(hilbert(DummyMat))));
    h_3(2)=plot(x, cat(2, flipdim(P_Filt(i, 2:end), 2), P_Filt(i, :)), 'r-', 'LineWidth', 1);
  end
  set(gca, 'YLim', [-60 0] + ceil(max(P_Filt(1, :))/5)*5)
else  
  h_3(2)=plot(x, 20*log10(max(abs(hilbert(cat(2, flipdim(P(:, 2:end), 2), P)*p_0/(1 + Sigma))))), 'r-', 'LineWidth', 1);
end
box on
xlabel('\bf{r} \rm{[mm]}')
ylabel('\bf{p} \rm{[dBPa]}')
title('maximum amplitude of the envelope')
